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I Introduction 

Self-assembly processes in patchy colloids represent one of the most striking examples where experimental methodologies and 
theoretical tools have progressed in parallel within a relatively short time scale^l While the former have been addressed 
elsewhere in this volume and in recent reviews 45 , in this contribution we will address the latter and, more specifically, the main 
methodologies that have been envisaged over the years to tackle the computation of the phase diagrams and phase transitions 
from one phase to another in dispersions of new-generation colloids, i.e., particles interacting via non-spherical potentials. 
Indeed, chemists and material scientists are starting to gain control on the shapes 6 and local properties of colloids. Hard cubes, 
tetrahedra, cones, rods as well as composed shapes of nano or microscopic size have made their appearance in the laboratories, 
and are becoming available in bulk quantities. Patterning of the surface properties of these particles 7 9 provides additional 
degrees of freedom to be exploited by scientists to engineer materials with peculiar properties. Patches on the particle surface 
can be functionalized with specific molecules 10 1 1 (including DNA single strandsS3Sl) to create hydrophobic or hydrophilic areas, 
providing specificity to the particle-particle interaction^^. 

Statistical physics provides a very rich and flexible toolbox to study thermo-physical and structural properties of complex 
fluids^, especially when coupled with the most recent and powerful computing techniques devised to deal with systems 
with a large number of degrees of freedorrt^USI While simple liquids and conventional colloidal systems have a long and venera- 
ble tradition^, theoretical studies on patchy colloids are relatively new, as in the past it was always tacitly assumed that even the 
unavoidable inhomogeneities in their surface composition could be neglected at a sufficiently coarse-grained scale. This is not 
the case, however, for patchy colloids that have surface patterns, chemical compositions and functionalities that are explicitly 
meant to be inhomogeneouJ^IHl Hence, the corresponding pair potentials describing inter-particle interactions depend on their 
relative orientations besides distances, and the analysis clearly becomes more complex. This is, however, by no means an insur- 
mountable difficulty, as several analytical and computational techniques have been devised in statistical physics to cope with the 
orientational dependence of the potential d 15 ! 16 ! 

In this Chapter, we shall discuss some of them in the framework of a particular pair potential that can be reckoned as a reasonable 
compromise between the complexity of the real interactions, and the necessary simplicity required to keep the analysis amenable. 
The basic idea of the model is built on the hard-sphere model, by providing a fraction of the surface sphere with a square-well 
character. This attractive region can be either condensed into a single large patch, or distributed over two (or more) patches 
symmetrically placed over the surface. Different spheres then interact via a square-well or hard-sphere potentials depending on 
their relative orientations and distances. 

This model was proposed in 2003 by Kern and FrenkeP and has ever since attracted considerable attention. There are two 
main reasons for this. On the one hand, the model is very flexible, as both the size and the number of the patches can be 
independently tuned, and this allows to mimic several different physical situations ranging from nanocolloids with more isolated 
attractive spots 21 to globular proteins with large regions of solvophobic exposed surfaced^!! On the other hand, the phase 
diagram obtained from the model can be directly compared with those obtained from experiments, as recently shown in several 
case^H 2 -^. In addition, the model displays some unusual features that can be paradigmatic for more complex systemsSSEl 
The aim of this Chapter is to introduce the main theoretical techniques to the evaluation of the phase diagram. This includes 
various Monte Carlo techniques (Section IV i, integral equations (Section|V]i, and perturbation theory (Section VI I. The level is 
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FIG. 1 : Sketch of one-patch parchy particles as modelled by the Kem-Frenkel potential. The surface of each sphere is partitioned into an 
attractive part (color code: blue) and a repulsive part (color code: white). The unit vectors fti and n 2 define the orientations of the pathces, 
whereas the vector f 12 joins the centers of the two spheres, from sphere 1 to sphere 2. The particular case shown corresponds to a 50% fraction 
of attractive surface (coverage \ — 0.5). 

intended to be pedagogic, with the main ideas behind each method outlined for non-experts in the field. Emphasis will be given 
on the calculations of thermodynamic quantities necessary for the phase diagram analysis, and hence a number of additional 
important results related to structural properties and other thermodynamical probes will be omitted. 

II The Kern-Frenkel model 

Consider a set of N identical hard-spheres of diameter a in a volume V at temperature T suspended in a microscopic fluid. When 
the surface of the spheres are uniform with no other interactions as their steric hindrance, the model has been often employed as 
a paradigm of sterically stabilized colloidal suspensions in the limit of high temperature or good solvent. 

As discussed elsewhere in this volume, colloids that are envisaged as elementary building blocks for the self-assembly process, 
are patchy colloids with different philicities in different part of the surface 45 . This means, for instance, that one fraction of the 
surface may be solvophilic and the other solvophobic. In solution, the solvophobic part will tend to avoid contact with the solvent 
and hence will act as an effective attractive force in the presence of another solvophobic patch laying on a different sphere. 
One can then consider the following model that was introduced in 2003 by Kern and Frenkel in the present forrrP^I but it is worth 
remarking that the idea of considering hard spheres and decorating them with patches of various forms and patterns dates back 
to much earlier, and several earlier versions in different fields can be considered as its ancestors^H^H. 

A circular patch is attached to the surface of each sphere, as depicted in Figure[T] with the central position of the patch identified 
by the unit vector n, and its amplitude measured by the angle 9q. Unlike the case of uniform sphere the interactions among 
spheres is anisotropic as it depends on the relative orientation of the unit vectors on each spheres with the direction connecting 
their centers. Then the idea is that two spheres attract each other if they are within the range of the attractive potential, with the 
corresponding attractive patches on each sphere properly aligned. 

If ni and n2 are the unit vectors associated with each patch on spheres 1 and 2, and f 12 is the direction connecting the centers of 
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the two spheres, then the interparticle potential reads 

$(12) = 0o (r 12 ) + (12) , (1) 
where the first term is the hard-sphere contribution 




0o (r) = { ■ ( 2) 

and the second term 

$i(ni,n 2 ,ri 2 ) = S w (n.a) * ("l, n 2 ,r 12 ) (3) 
is the orientation-dependent attractive part which can be factorized into an isotropic square-well tail 

!— e, g < r < Act 
(4) 
0, Act < r 

multiplied by an angular dependent factor 

I 1, if ni • fi2 > cos O and -n 2 • f 12 > cos 9 
*(ni,n 2 ,ri 2 ) = < . (5) 

I 0, otherwise 

Here a is the sphere diameter, (A — l)er is the width of the square-well interaction and e its depth. 28q defines the angular 
amplitude of the patch. The unit vectors hi(oJi), (i — 1, 2), are defined by the spherical angles to* = (8i, </?,-) in an arbitrarily 
oriented coordinate frame and ?i 2 (0) is identified by the spherical angle il in the same frame. Reduced units, for temperature 
T* = ksT/e, pressure P* = f3P/p and density p* = per 3 , will be used throughout, with ks being the Boltzmann constant. For 
future reference, we also introduce the packing fraction rj = np* /6. 

The coverage x is me fraction of attractive surface on the particle. \ can t> e related to the patch half-width 9 as 

x = h(n u n 2 M 1/2 =1^° (6) 

\ V / / UJ1UJ2 

where we have introduced the angular average 

(-..). = il^-. (7) 



III The tools of statistical physics 

Statistical physics has developed a number of different theoretical approaches to compute the thermophysical properties of a 
fluidCK] In order to compare with experiments, we are most interested in the computation of the fluid-fluid and fluid-solid 
phase diagram on the one hand, and on the specific mechanism driving aggregation, and hence self-assembly, on the other 
hand. Among this arsenal of different available techniques, here we will review three different methodologies that were recently 
exploited in the framework of the Kern-Frenkel model. These are Monte Carlo simulation ;) 17 ! 18 !, integral equation theories^ and 
thermodynamical perturbation theories^. 
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Monte Carlo simulations are undoubtably one the most efficient ways to accurately compute the properties of a model fluid. 
As discussed in more details below, the main limitations of simulations are that they can be very demanding from a computa- 
tional point of view, especially for sufficiently realistic potentials, and that they are unable to distinguish metastable from stable 
equilibrium states. On the other hand, they provide virtually exact estimates of all static quantities of interest. Several improve- 
ments have been proposed over the years, some of them triggered by the problems discussed in this contribution, so that the 
methodology is very well established and by now extensively reviewed and described in detail in several books (see Refs. 1171181 
and references therein). The case of patchy colloids, however, is relatively recent, although it builds upon previous established 
procedures on other complex fluids. 

Integral equation and thermodynamics perturbation theories are two of the main methodologies from the toolbox of Statistical 
Physics that are at the basis of our current understaing of simple and molecular fluid d 15 l 16 4 In spite of their known drawbacks and 
shortcomings, they are known to provide reliable predictions for both structural and thermodynamical properties each in their 
own domain of applications. Their applications to patchy colloids is a rather natural, albeit not trivial, extension of formalisms 
already developed in the last two decades for molecular fluidP. As it will become clear, they both become particularly attractive 
in view of the large computational effort involved in Monte Carlo simulations. In addition, they are able to access to some details 
and nuances that are not easily accessible by other methods. 

IV Monte Carlo simulations 

The aim of Monte Carlo simulations is the computation of thermodynamic quantities by performing an average over a suitable 
ensemble of microstates. The choice of the ensemble is dictated both by the quantities to be computed and by the specific system 
under investigation, for which one ensemble can be more convenient than the others. Below we review the most interesting 
techniques that have been used to calculate phase diagram of patchy colloid models. 

A Canonical NVT and NPT methods 

Simulations in the NVT (isothermal-isochoric) and NPT (isothermal-isobaric) ensembles are probably the most common 
example. In these ensembles, the number of particles N, the temperature T and the volume V (NVT) or the pressure P (NPT) 
are held constant. The Markov chain in configuration space is constructed via a sequence of translational and rotational moves, 
accepted with an appropriate probability that depends upon the change in potential energy and T. In the NPT case, the volume 
is also varied. With a proper choice of the acceptance probability, the system first evolves toward equilibrium and then starts to 
sample equilibrium configurations with the Boltzmann statistical weight. The equilibration process can be rather long, especially 
in cases where kinetic traps are present (as in the vicinity of gel or glass transitions) or when an activation barrier needs to be 
overcome. This last case arises when the system is metastable with respect to a lower energy phase or when it organizes into 
mesoscopic structures and specific self-assembly processes involving large numbers of particles are requested. The approach to 
equilibrium can be monitored by focusing on the time evolution of collective properties (e.g. the potential energy, the density, 
the pressure). Since equilibration can be rather slow, it is highly recommended to make use of a logarithmic time scale when 
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searching for a drift in the time dependence of the investigated property. 

When a sufficiently large number of statistically independent equilibrium configurations have been generated and stored, all 
possible structural (static) information can be calculated. Typical quantities that are computed are the total energy U, the radial 
distribution function g(r) and the structure factor. In the case of anisotropic systems, such as patchy particles, the orientational 
dependence of the structural properties needs to be evaluated as well. The center-center g(r) is indeed not sufficient to evaluate 
the average potential energy or the pressure, at variance of the isotropic case where calculation of U or P from the g(r) usually 
simply requires a one-dimensional integration. 

In the case of hard bodies or in the presence of step-wise potentials (e.g. the square-well potential), direct evaluation of P in 
NVT simulations is in principle possible^^but not straightforward. To evaluate the equation of state, i.e., the relation between 
density and P at constant T, the NPT ensemble is often preferred in this case. 

Various additional improvements can be (and are) used to improve the convergence of the scheme in a way that will be described 
in each specific example. 

B Gibbs ensemble method 

A convenient scheme was devised by Panagiotopulos 38 to specifically address the problem of a direct evaluation of the gas- 
liquid phase coexistence by Monte Carlo simulations. This is known as the Gibbs Ensemble Monte Carlo (GEMC) method. N 
particles are partitioned into two distinct simulation boxes. In addition to intra-box translational and rotational moves, particle 
and volume swap moves (keeping both the total number of particles and the total volume fixed) are proposed and accepted with 
the appropriate probability 11 . In this way the two coexisting phases are simulated without the intervention of a interface between 
them. When convergence is reached, the densities in the two boxes provide the value of the coexisting densities of the liquid and 
gas phase. It should be pointed out that the GEMC method becomes inefficient when the density of the liquid phase becomes 
large, since the probability of inserting a particle with in a favourable state — i.e., not overlapping with any other — becomes 
extremely small. 

For the specific case of KF particles discussed later in this Chapter, GEMC simulations have been performed for a system of 
1200 particles in a total volume of (16er) 3 . On average, the code attempts one volume change every five particle-swap moves 
and 500 displacement moves. Each displacement move is composed of a simultaneous random translation of the particle center 
(uniformly distributed between ±0.05cr) and a rotation (with an angle uniformly distributed between ±0.1 radians) around a 
random axis. 

C Grand-canonical ensemble \iVT 

In the neighborhood of the gas-liquid critical point the free-energy barrier separating the two phases becomes comparable to 
the amplitude of the thermal fluctuations of of the relatively small systems that can be accessed in simulation. In this case, the 
GEMC method cannot be used to investigate the system since it becomes size effects and spontaneous swapping of the average 
densities between the two boxes. 
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A precise evaluation of the critical parameters (density and temperature) can be obtained performing simulations in the grand- 
canonical [iVT ensemble 17 , where density fluctuations are accounted for at fixed volumes and temperature, coupled with the 
finite-size scaling analysis envisaged by Wilding 39 . MC simulations in the gran canonical ensemble are implemented by per- 
forming trial insertions and deletions of particles, besides trial displacements and rotations. The critical parameters of the system 
can be extracted by matching the calculated distribution of density fluctuations to the expected distribution at the critical point, 
a feature which is largely system independent. 

In the implementation of the grand-canonical simulations to the KF model reported later on, one insertion/deletion attempt was 
performed, on average, every 500 trial translational/rotational displacements. 



D Fluid-Solid coexistence: thermodynamic integration 

To compute numerically the free energies of the fluid and the crystals and their coexistence lines it is possible to resort to 

thermodynamic integration methods. Details of this procedure were recently given in a detailed review by Vega et aP3. 

The starting point of the procedure requires the identification of a state point in the pressure-temperature plane where two phases, 

I and II, share the same chemical potential, i.e., T) = /j,u(P,T) . The chemical potential of the fluid can be computed 

by thermodynamic integration using the ideal gas as a reference state, and by integrating the equation of state, P(p), at fixed 

temperature 



N °* > J„ p> 

where F/N is the Helmholtz energy per particle. The first term on the right-hand-side of Eq.([8]l is the ideal gas part and depends 
upon the system dimensionality. The chemical potential can then be recovered as 

Plx{P{p),T) = PF(P ( phT) +0P( P )/p. (9) 
To compute the chemical potential of a crystal one can perform thermodynamic integration at fixed density and T using an 



ideal Einstein crystal as the reference system. This method, known as Frenkel-Ladd procedur e! 17 ! 40 ! is very efficient and by now 
standard. Integration of the crystal equation of state provides a way to evaluate the chemical potential at different T and P. The 
pressure at which the chemical potential of the fluid and of the crystal are identical along an isotherm provides the coexisting 
pressure at the selected T. 

Starting from a coexistence point, coexistence lines can finally be inferred by using Gibbs-Duhem integration, as described by 
KofkeP, numerically integrating over the Clausius-Clapeyron equation. 



V Integral equation theories 
A General scheme 



At first, let us consider simple fluids first where the particles can be regarded as spherically symmetric. All thermodynamic 
properties of such fluids can be straightforwardly computed from the radial distribution function g(r). In integral equation 
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theory 15 , the strategy to infer the thermophysical properties of a fluid hinges on the calculation of the total correlation function 
h(r) = g(r) — 1 that, in turn, is related to the direct correlation function c(r) by the Ornstein-Zernike (OZ) equation 



h(r) =c(r) + p J dr'c(r') h(\i 



(10) 



Once that h(r) is known, all statistical and thermodynamical properties can in principle be computed. However, as h(r) de- 
pends upon the unknown quantity c(r), an additional equation involving both quantities is required for the solution. Unlike 



equation (lOi, which is exact, the second equation always involves some approximation. This gives rise to some well known 
thermodynamic inconsistencies, that are the main shortcomings of this method, and that may severely limit its applicability. 
The second relation between h(r) and c(r) also involves the two-body potential <fi(r), and can be cast in the general form 

c(r) = exp [-P<p (r) + 7 (r) + B (r)] - 1 - 7 (r) (11) 

where 7 (r) = h (r) — c (r) is an auxiliary function. Although this equation is again exact in principle, it involves the bridge 
function B(r) that in general depends upon higher body correlation functions, so in practice an approximation (closure) is 
always necessary. The quality of the results obtained will depend crucially on the reliability of the involved approximations; 
several closures have been proposed over the years with their pros and cons well classified and under control. Among them, the 
reference hyper-netted chain (RHNC) stands out as an optimal trade-off between simplicity and precision of its predictions, and 
this is the one that will be the object of the present Chapter. 

Having closed the systems of two equations in two unknowns (h(r), and c(r)), the system may then be solved iteratively with 



the convolution appearing in the OZ equation ( 10 1 simplified in Fourier space, as 



h(k) = (12) 
1 — pc (k) 

h(k) and c{k) being the Fourier transform of h(r) and c(r) respectively. 

The RHNC closure was introduced by Lade 42 for spherical potentials and later extended to molecular fluid d 43 ! 44 ! In the RHNC 



closure, one replaces the exact B{r) appearing in ( 11 1 by its hard-sphere counterpart Bo(r), that is the only system for which 
a reliable expression (the Verlet- Weiss expression 4 -^ is available. Rosenfeld and AshcrofP demonstrated that the effectiveness 
of the reference system could be magnified by treating its parameters as variables to be optimized in some fashion. It is in fact 
possible to determine them via a variational free energy principle 47 that enhances internal thermodynamic consistency. With the 
effective hard sphere diameter <jq suitably chosen in this way, the RHNC has been shown to provide rather precise estimates 
of the chemical potential and pressure, that is the two crucial thermodynamical quantities needed for the calculation of phase 
diagrams. 

The case of anisotropic potentials is significantly more complex from the algorithmic point of view, but the philosophy behind 
the methodology is identical. The procedure hinges on a remarkable piece of work carried out by Fred Lado in a series of 
p a p er ^43|44|47]j n me framework G f molecular fluids and more recently adapted to the case of patchy colloids. Here we just sketch 
the idea, referring to Refs. |48Tl50l for details. 



The angular dependent counterparts of Eqs.( 10 1 and ( 1 1 1 are given in terms of 7(12) = h(12) — c(12), and are 



7(12) = £ y , dr 3 d W3 [7(13) + c(13)]c(32), (13) 
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/ \ Fourier transfoiTn ~ / 7 \ 

c(r) y c(k) 



Closure 



OZ equation 



/ \ Inverse Fourier transfoiTn „ / , \ 

7(r) < 7(fc) 



TABLE I: Schematic flow-chart for the solution of the OZ equation for isotropic potentials, 
for the OZ equation, and 

c(12) = exp [-/?$> (12) + 7 (12) + B (12)] - 1 - 7 (12) . (14) 

for the closure equation. Again, the RHNC approximation amounts to assume B{V2) = B (ri 2 ), but clearly this is a much 
more drastic approximation in the present case, as the real 5(12) depends on the relative orientations of the particles whereas 
the reference B$(r 12) does not. As a result, one might expect the results to be less precise in this case as compared with the 
isotropic fluid counterpart. 



B Iterative procedure 

As discussed before, the iterative procedure in the case of the spherical isotropic potential is very simple and outlined in Table 
I. It requires a series of transformation to and from Fourier space, where the solution of the OZ equation is more conveniently 
carried out in view of Eq.(|T2]> The iterative solution of the angular dependent Ornstein-Zernike (OZ) equation ( p"3j ) along with the 



approximate closure equation ( 14 1 again requires a series of direct and inverse Fourier transforms between real and momentum 
space involving the bridge function _B(12), the direct correlation function c(12), the pair distribution function <?(12), the total 
correlation functions h(12) and the auxiliary function 7(12) = h(12) — c(12). 

In addition to this, however, the orientational degrees of freedom introduce additional direct and inverse Clebsch-Gordan (CG) 
transformation between the coefficients of the angular expansions in different frames^. Two important examples are the so- 
called "axial" or "molecular" frames, with z = r.y in real space, and the {k} representation with z = k in momentum space, 
because in those representations some of the calculations become particularly simple. This set of transformations also allows 
the definition of the correlation functions (in particular the g(12)) in an arbitrarily-oriented axes. We further note that, in the 
presence of an anisotropic potential, the Fourier transform is implemented through a Hankel transform, and that the cylindrical 
symmetry of the angular dependence included in the Kern-Frenkel potential of Section|II](when the number of patches present 
on each sphere is one or two) allows us to use the simpler version of the procedure for linear molecules. 

All necessary equations can be found in Ref. 16 that is the standard reference for this topic, and only the most important ones 
will be reported in the following. 

The resulting scheme is illustrated in Table II and is the extension of the isotropic case given in Table I 49 . Consider the expansion 
in spherical harmonics of the auxiliary function 7(12) in the axial frame 

7( 12 ) = 7 ( r , Ui, ^2) = 4tt ^2 Jhhm (r) Y hm (uii) Yi 2fh (uj 2 ) , (15) 

li ,12 ,tti 
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/ t 7 7 \ Hankel transform _, / , , , , <. Inverse CG transform ~ / , \ 

c(r;hl 2 l) > c(k;hl 2 l) > c hhm {k) 



CG transform 

Chl 2 m if) 

Expansion inverse 

Closure 

7(r, 

Expansion 



[ihhm (r)] 



old 



Iterate 



OZ equation 

Ihhm (k) 

CG transform 

Inverse Hankel transform 

l{r;lxhl) 

Inverse CG transform 

[ihhm (r)] 

new 



TABLE II: Schematic flow-chart for the solution of the OZ equation for the Kern-Frenkel angle-dependent potential. See Section [V] for a 
description of the scheme 



where m = —to, and its inverse reads 

Ihhm (r) 



in 



dujxduj^l (r, Wi, w 2 ) Y t * m (wi) Y,* rn (w 2 ) . 



(16) 



Eq.( 16 1 provides the initial set [7; 1 i 2m ('")]oid described as the initial point in Table II, whereas Eq.( 15 i yields the next term in the 
iteration map j(ri2, Wj., £^2). 

Using the aforementioned RHNC closure approximation -6(12) = B (ri 2 ), the bridge function is constructed and then inserted 



in Eq.( 14 1 to get c(r, ui, Eq.( 16 1 with the replacement 7 — > c is then exploited to infer the corresponding axial coefficients 



Qi/ 2 m('"l2) 

The next step is to implement a Clebsh-Gordan transform in direct space, in order to transform from axial coefficients c; 1 / 2m (r) 
where z = f to space coefficients c (r; M2O in an arbitrarly oriented frame. The necessary expressions are the equation pairs 



c(r;hl 2 l) 



-in 



21 



1/2 



-j C (hhh mmO) c h i 2m (r) 



where the C (I1I2I', mmO) are Clebsch-Gordan coefficients, with the inverse transform given by 



Q 



ihm(r) = y^C(lil 2 l;mmO) 



21 + 1 
in 



1/2 



c(r;hl 2 l) ■ 



(17) 



(18) 



and the coefficients c(r; ii^O are then given by Eq.( 18 1 



The last required tool is the Fourier transform of the radial parts, that is a Henkel transform as given by the pairs 



POO 

c{k;lxhl) = 47ti' / dr r 2 c(r;lil 2 l) ji (At) 
^0 



with the inverse transform reading 



c(r;hl 2 l) 



2n 2 \ l 



dk k 2 c (k; hhl) ji (kr) , 



(19) 



(20) 
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where ji(x) is the spherical Bessell function of order I. By using "raising" and "lowering" operations on the integrands these 
are finally evaluated with jo(kr) = s'mkr/kr kernels, for I even, and j_i(kr) = coskr/kr kernels, for I odd. Fast Fourier 
Transforms are programmed for both instances. (Note that Hankel transforms X (k; l\l 2 l) are imaginary for I odd while for I 
even they are real.) 



Having obtained c(fc; l\l 2 l) from Eq.( 19 1, we have then reached the turning point in Table II, from which the returning part 
can then be started with a parallel sequence of operation in Fourier space. These include, a backward Clebsch-Gordan trans- 
formation to return to an axial frame in k space and get 2/ 1 ; 2l „(fc); a OZ equation in k-space to get 7; 1 ; 2m (fc), followed by a 
forward Clebsch-Gordan transformation, and an inverse Hankel transform, to find j(r; hl 2 l). A final backward Clebsch-Gordan 
transformation, brings a new estimate of the original coefficients [7i 1 j 2m (r)] new . This cycle is repeated until self-consistency 
between input and output 7/ 1 z, m (V) is achieved as before. 

Note that the OZ equation required in {k} representation's expressed in terms of the axial expansion coefficients of the trans- 
formed pair functions, is given by the following matrix form 



(k) = (-1)">^ [%l 3 m 0) + C h i 3m (k)] Q 3 ; 2m (fc) , 
i 3 =m 



(21) 



C Thermodynamics 

Once that the correlation function h(12) (and hence the distribution function ,g(12) = h(12) + 1) is known, the excess free 
energy can be computed a d 44 ! 49 ! 



TV 



N 



PF 2 
N 



PF 3 
N ' 



where 



PF X 


1 


N 


= ~2 P 


(3F 2 


1 


N 






PF° 


N 


N 



dr 12 ( -h 2 (12) + h (12) - g (12) In \g (12) e mi2) 



(27T) 



In Det 



I + (-l) m Ph m (k)] - (-l) m pTr [h m (*)] } , 



\p J dr 12 {[g(12)-g (12)]B (12)) UiU 



(22) 

(23) 
(24) 

(25) 



In Eq. (24 1, h m (fc) is a Hermitian matrix with elements /i; 1 / 2m (fc), l\, l 2 > m, and I is the unit matrix. The last equation, for 



F 3 , directly expresses the RHNC approximation. Here is the reference system contribution, computed from the known free 
energy F<? x of the reference system as F§ = F c ° x — F® — F 2 , with F® and F 2 calculated as above but with reference system 
quantities. 



The bridge function _Bq(12) appearing in (25 i is the key approximation in the RHNC scheme, since it replaces the unknown 



bridge function 5(12) in the general closure equation ( 14 1. This is taken from the Verlet-Weiss expression of the hard-sphere 
model 45 , as anticipated, in view of its simplicity and of the fact that it works reasonably well for the case of the square-well as 
we will see, but with a renormalized diameter ctq for the hard-sphere that is selected by enforcing the variational condition 47 



7 r / \ / \i dB HS (r;cr ) 
(> I dr [gooo (r) - g m [r; a a )\ — = 0. 



(26) 



12 

From the free energy F, one of course compute all thermodynamics following standard procedures. For the computation of the 
phase diagram, the pressure and chemical potential are required at any fixed temperature. 
The virial pressure P is obtained as^ 

1 / N N B \ 1 f I 8 \ 

P = pkBT -W\^^ rij d^* iij) ) =P k * T -- & P l ] rfr 12 ^. 9 (12)r 12 — $(12)^^ . (27) 



that in turns can be cast in the form involving the cavity function 2/(12) = g(12)e ,3 *( 12 ) with the result 49 

l + lnp^llyia^Lo^e^^A - A 3 (y (Xa, Wl , a*) [e**^'*' ^ - ll ) j, (28) 



£P 1 , 2 
P 

As already remarked, one of the main advantages of the RHNC closure stems from the fact that the calculation of the chemical 
potential does not introduce any approximation in addition to that already included in the closure. It can be obtained from the 
thermodynamic relation 

N p 

that was already used in Eq.Q. 

Finally, we note that the ideal quantities for the free energy, the virial pressure and the chemical potential are 

/3^ id ,_, .3, „ l3P id 



N 



m (pA 3 )-l ^ = 1 /3 Mid = m ( p A 3 ) (30) 



VI Barker-Henderson perturbation theory 

Another powerful method to access thermophysical properties of a fluid is thermodynamical perturbation theory, that directly 
extracts the free energy of the system from the knowledge of the free energy F of a reference fluid (the hard-sphere fluid in 
the present case). This is a well-known techniques in several fields of physics, including simple 15 and complex 16 fluids, and 
hinges on the fact that often the system under investigation is not very different from the reference one, so that an expansion in 
this perturbation term is a reasonable approximation. Under these conditions, the results are expected to be rather reliable, even 
when stopping the expansion at the lowest orders. 

In the square-well fluid case the analysis has been carried out in details in the late starting from the pionieering work by 

Zwanzig 53 , and recently extended to the patchy cas d 23 l 54 l 

Working in the grand-canonical ensemble as this is the most convenient one^l, we assume the total potential U to have the 
following form 

U 7 (l,...,N) = U (l,...,N)+ 1 U I (l,...,N) (31) 

i<j i<j i<j 

where f7o(l, ■ ■ • , N) — J2i j ®o(ij) is the unperturbed part and . . . , N) = J^i j i s tne perturbation part. Here 

< 7 < 1 is used as perturbative parameter. Note that when each coordinate i includes both the coordinate and patch 
orientation n,, so that i = (r^ , n.j), then the expression is valid also for the Kern-Frenkel model 23 54 . For simple fluids, then 
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Yi only. Introducing the following short-hand notation 



l,....N 



(...) 



N 



(32) 



for the integration over all particle coordinates, the grand-canonical partition function 

+00 

2, - E 



+°° e PnN 



N\A™ J, N 



(33) 



N=0 ■' T 

(here At is the de Broglie thermal wavelength, and fi 7 is the grand-potential) can then be used to obtain an expansion of the 
Helmholtz free energySl. 

'dF, t 



F~ 



2\ 7 \ 9 7 2 



(34) 



7=0 



that is valid for arbitrary 7. 

Taking the derivative of In Q 7 at fixed chemical potential p, one has, using Eq.( 3 1 

d 



97 



In e_ 



1 



d 



2 J 1,2 d-f 



[-/3$ 7 (12)]p 7 (12) 



where 



+00 



py(l...h) 



-pu-, 



(35) 



(36) 



N=h v ~ ' T 

When 7 = 1, this yields the free energy correct to first order in the expansion Eq.([34|. 

The second order correction is far more laborious. Indeed, the extension of this analysis involves higher orders correlation 
functions and this forces additional approximations to come into pla y ! 51 ! 52 !, mus humpering its practical utility. In 1967, Barker- 
Henderson gave a much simpler recipe^ that was found to be rather effective in predicting the phase diagram of the square-well 
fluicP and was recently extended to the Kern-Frenkel case®! 

The final result for arbitrary angular dependent potential, correct to second-order, reads^ 



f3F_ = /?F f3F 1 f3F 2 
N N N N 



where 



and 



P N / drg (r) (^(r.fi.ui.oto))^ 



PF 2 



drg Q (r) ([/?$! ( r ,Sl, wi,wa)]' 



(37) 



(38) 



(39) 



In the particular case if the Kern-Frenkel potential, one obtains 23 54 



N 



12/7 

~^3~ 



drg (r) (09 (12)) 



(40) 
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and 



N 



6r/ f dr\ 



a 3 \dPJ 



dvg (r)4 w (r) (1^(12) 



wl,W2 



(41) 



Here Pq = f3Po / p is the reduced pressure of the HS reference system and go ( r ) the corresponding radial distribution function. 
As before, the pressure and chemical potential can be derived from the reduced free energy per particle (3F/N, using the exact 
thermodynamic identities 



P 

pp = 



n 



d fpF 



dr) 



N 
p_F 
N 



(42) 
(43) 



The phase diagram in the temperature-density plane can then be computed using a numerical procedure that will discussed in 



connections with Eqs.(44i and (45 i below, as it is identical to the one used for integral equations. 



VII Calculation of the fluid-fluid coexistence curves for integral equation and pertubation theory 

The common feature of integral equations and thermodynamic perturbation theory is that one is able to obtain approximate 
expressions for both the pressure P and the chemical potential p as a function of the temperature T and the density p. In the 
presence of a fluid-fluid (gas-liquid) transition, both P and p will have well defined dependence on T and p in the gas and liquid 
branches, but not in the coexistence regions. Hence, in order to obtain the coexistence curves, and hence the phase diagram in 
the temperature-density plane two procedures are possible. First a graphical procedure where one plots the chemical potential 
versus the pressure at a given fixed temperature, and seeks the intersections between the gas and the liquid branches. An example 
of this procedure in the case of the square-well potential was given in Ref. 48 This procedure, however, is neither very practical 
nor very precise, as it involves qualitative deduction of the crossing points. Yet it can be used as a first preliminary estimate of a 
more precise numerical calculations as follows. 

For a fixed temperature T, one can then compute the pressure of the gas (colloidal poor) phase P g and of the liquid (colloidal 
rich) Pi phases, and the corresponding chemical potentials p g and pi. The fluid-fluid (gas-liquid) coexistence line then follows 
from a numerical solution of a system of non-linear equations 

P g {T, Pg ) = Pi{T, Pl ) (44) 
Ms (T,p g ) = pi (T,pi) (45) 

whose solutions are the gas p g and liquid pi densities associated with the coexistence lines. By plotting the resulting p g and pi as a 
function of T the coexistence curve can be constructed in the region where the transition occurs. This procedure will be followed 



in the analysis of the Kern-Frenkel fluid-fluid phase diagram for both integral equations (see Section VIII A i and thermodynamic 



perturbation theory (see Section VIII D i, but can also be exploited for the computation of the fluid-solid transition, as explained 
in Section lyTTTEl 
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FIG. 2: Fluid-fluid coexistence curves ksT/e versus pa for the SW fluid of range A = 1.5. Data from RHNC integral equation (squares) are 
contrasted with Monte Carlo simulations by Vega et alp^ (circles), and del Rio et alP^ (crosses). Adapted from Refl48l 

VIII Results 

A Fluid-Fluid coexistence curves from RHNC integral equation 



We start by reviewing the quality of the RHNC integral equation theory for the simple case of the isotropic square-well case, a 
model which has often be used as an effective potential (in the implicit solvent representation) to model colloidal particles inter- 
acting via depletion interactions^^!. Fig j2] shows the RHNC predictions for the gas-liquid coexistence from Ref 48 contrasted 
with Monte Carlo simulations on the same system carried out by different groupJ^H. 

Although the integral equation results appear to reproduce reasonably well those from numerical simulations, two features are 
worth noticing. 

First, thermodynamic inconsistencies associated with the approximate nature of the closure are at the origin of the incorrect 
evaluation of the pressure and chemical potentials and hence of the exact location of the coexistence lines. This is an unavoidable 
feature of all integral equations, and its origin and effects are well known. For most of the cases, in fact, the whole critical region 
is inaccessible to integral equation theories. 

A second additional point is related to the non-linear nature of the self-consistence procedure, and gives rise to numerical 
instabilities that may or may not be controlled depending on the considered state point. As a general rule, lower temperatures 



and higher densities are harder to converge, and hence for some points the solution of the system of Eqs. (44 1 and (45 1 might 
not even exist. 

These considerations are even more compelling in the more complex case of the Kern-Frenkel fluid where condensation takes 
place at lower T associated with the involved lower coverages. Hence one might expect an agreement with respect to numerical 
simulations not better than what has been reported for the isotropic square-well case. This is indeed the case, as shown in Fig(3] 
for the representative example of single-patch particles with coverage \ = 0.8, for which 80% of the surface has a SW character, 
condensed into a unique patch. The width of the square-well has been selected to be A = 1.5 as before, so that the limit of full 
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FIG. 3: Fluid-fluid coexistence curves ksT/e versus per 3 for the Kern-Frenkel fluid with coverage \ = 0.8 and range A = 1.5. Data from 
RHNC integral equation (circles) are compared against MC simulations (squares). Adapted from Refl49l 

coverage gives back the result of Fig j2] and the value \ = 0.8 has been selected to be half-way between the fully occupied fluid 
and half-coverage that has peculiar behavior, as will be discussed shortly. 

In Fig[3]we also report the coexistence lines and critical points using Gibbs ensemble and Gran Canonical Monte Carlo simu- 
lations, as outlined in Sections [IV B | and |I V C | The former were used to evaluate coexistence in the region where the gas-liquid 
free-energy barrier is sufficiently high to avoid crossing between the two phases, whereas the latter was used to locate exactly 
the critical point. The details of the analysis can be found in Refs. [48](and in Ref.|50]for the corresponding two-patch case). 
The comparison between RHNC integral equation and MC simulations for the Kern-Frenkel model with intermediate coverage 
X = 0.8 displayed in Fig|3]shows that indeed the agreement is comparable with the square-well case as anticipated. 
The Kern-Frenkel model offers the possibility to continuously change the coverage interpolating from the isotropic square-well 
to the symmetric Janus-like potential, when the coverage moves from \ — 1 to X — 0.5 (see Fig. |4]-(a)). To investigate how 
the phase diagram of Janus particles arises, we calculate how the gas-liquid coexistence is modified on progressively reducing 
X- Fig. |4]-(b) shows MC simulations (Gibbs ensemble for the coexistence lines and gran-canonical fiVT for the critical point) 
results for the gas-liquid phase coexistence for several \ values, extending the original data by Kern and FrenkeP. A progressive 
shrinking of the coexistence region to lower temperatures and densities accompanies the decrease of the coverage. Consistently, 
as the coverage decreases, both the critical temperature and critical density decrease (see Fig|7]). This is not surprising, in 
view of the fact that the coverage is a measure of the attractive interactions intensity (by controlling the maximum number of 
nearest neighbors contacts^®!) that dictates the value of T c and (perhaps more indirectly) p c . It can be explicitly shown that this 
shrinking of the coexistence region is a non-trivial one, and that it cannot be inferred by a simple scaling of either the temperature 
or the density. Up to x = 0.6 coverage, however, the morphology of the curve appears to be the standard one, with the gas and 
liquid coexistence lines widening on cooling down. 

The half-coverage x = 0.5 is known as the Janus limit and will be discussed in the next Section, as it displays an interesting 
unconventional behavioi^. 



17 



(a) Square Well 



Janus 



©ooo 



*=1 



Z = 0.8 



^ = 0.6 




FIG. 4: (a) Cartoon of a one-patch particle with the coverage (fraction of attractive surface, depicted in blue) changes from one (square well) to 
one half (Janus), (b) Fluid-fluid coexistence curves fcaT/e versus pa 3 for the Kern-Frenkel fluid with descreasing coverages from x = 1-0 (the 
SW case) to the Janus limit \ = 0.5. The range is always set to A = 1.5. Data are from Gibbs ensemble MC simulations for the coexistence 
lines and from grand-canonical MC simulations for the critical points. This is the one-patch result: for the two-patches counterpart, see Fig|6] 
Adapted from Refl3T1 



B The Janus limit 



The value of x = 0.5 plays a special role in the Kern-Frenkel model with a single patch. This can be seen by following the 
change in the coexistence line locations as the coverage is reduced from the fully isotropic square- well fluid (x = 1.0) to the 
Janus case (x — 0.5), as discussed in the previous section. 

In Fig. [5] (a) the Janus case x = 0.5, already reported in Figure |4] is shown by magnifying the very narrow region where the 
transition occurs. 

The difference with the standard phase diagram appears rather clearly. As the fluid is cooled down to lower and lower tem- 
peratures, the coexistence region shrinks, contrary to the standard case, and the two coexistence lines appear to approach one 
another at sufficiently low temperatures. The origin of this anomaly has far reaching consequences that have been discussed 
m detailPEl. The crucial point is that at low temperature and density, monomers tend to aggregate in micelles and vesicles 
(see Fig[5]-(b)) bearing a well defined number of particles (10,41, . . .) so that all favorable contacts are saturated inside the 
aggregate ! 28 ^ 3 H This means that, in all cases, the clusters always expose the hard-sphere parts to the external fluid, and the 
condensation process is thus inhibited and eventually destroyed altogether. The specificity of these magic numbers can be also 
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Vesicle 

FIG. 5: (a) The anomalous phase diagram of the Janus limit (x = 0.5), adapted from Ref 28 (b) Cartoon of the aggregates which form in the 
gas phase on cooling, micelles and vesicles. The blue surface is attractive, modeled via a square-well potential. 

explained by means of a simple cluster theory 61 , and an even simpler approach 62 can be developed to mimic the onset of the 
re-entrant gas branch. Below the Janus limit, no evidence of fluid-fluid transition was found, in full agreement with the above 
interpretatiotPEl. 

The unconventional shape of the phase diagram arises from the particular energetic and entropic balance associated to the 
transition from the micelles gas to the liquid phase. It has been found that, differently from the usual gas-liquid behavior, the 
potential energy per particle is higher in the liquid phase than in the micellar gas phase, a consequence of the greater energetic 
stability of the micelles and vesicles as compared to the disordered liquid phase. Hence, despite the fact that the gas phase is 
stabilized by the translational entropy of the micelles, the coexisting liquid phase is more disordered than the gas one. Such 
unconventional entropic stabilization of the liquid phase arises from the orientational entropy, since particles are orientationally 
disordered in the liquid phase while they are properly oriented in the micelles gas phase. 
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FIG. 6: (a) Cartoon of two-patch particles, when the coverage (blue) varies from one (square well limit) to zero (hard sphere limit), (b) 
Fluid-fluid coexistence curves k^T/e versus pa 3 for the Kern-Frenkel fluid with two-patches and descreasing coverages from x = 1-0 (the 
SW case) to the value x = 0.3. The range is always set to A = 1.5. Data are from Gibbs ensemble MC simulations for the coexistence lines 
and from grand-canonical MC simulations for the critical points. This is the two-patches counterpart of Fig|4] Adapted from Refl50l 

C One versus two-patches 

It is interesting to investigate the effects of having the same attractive square-well coverage split in two parts, at the opposite poles 
of the sphere so that they are symmetrically distributed. This is the two-patche case, and again this case smoothly interpolates 
between the fully occupied isotropic SW case x = 1 an d the empty hard-sphere x = case. The corresponding fluid-fluid 
(gas-liquid) phase diagram is depicted in Fig j6] that should be contrasted with the single patch counterpart reported in Figj4] 
Two main differences are noteworthy. First, the case x = 0.5 does not appear to play any particular role, at variance with the 
single patch counterpart. Its coverage value corresponds to the so-called tri-block Janus case, and will be further discussed later 
orPSEI]. This can be attributed to the fact that at higher valences it is not possible to saturate all favorable contacts, and hence the 
condensation process always takes places in complete agreement with the previous interpretation of the Janus case. The second 
main difference is related with the fact that in the two-patch case coverages below x = 0.5 do exhibit fluid-fluid transition, as 
displayed in Figj6] where coverages as low as x = 0.3 are depicted. Below this value, the transition becomes metastable with 
respect to crystallization as explained in details in Ref.l50l 

The x dependence of the critical point can be inferred by plotting the reduced critical densities and temperatures, as a function 
of the different coverages, as shown in Fig|7j ! Clearly, the two-patch densities and temperatures lie always above the one-patch 
counterparts, thus supporting the interpretation that it is easier to form a fluid with higher valence at a given coverage, a feature 
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FIG. 8: The fluid-fluid phase diagram from BH thermodynamic perturbation theory (continuous lines) contrasted against numerical simulations 
(open symbols). Closed symbols refer to the BH critical points. Adapted from Refl54l 



that is related to the existence of the fluid-fluid transition even for low coverages. 



D Evaluation of the fluid-fluid coexistence curves from thermodynamic perturbation theory 

As in the case of integral equation theory, the fluid-fluid phase diagram for in the temperature-density plane can be computed 
even from Barker-Henderson (BH) thermodynamic perturbation theory^, as indicated in Section 



VII 



This is shown in Fig.[8]for the one-patch case, and compared with the same MC simulations used for comparison with integral 
equation theory. 

Given the simplicity of the theory, the accuracy of the BH theoretical prediction is rather remarkable. On the one hand, this 
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FIG. 9: The fluid-solid phase diagram from BH thermodynamic perturbation theory (continuos lines) contrasted against numerical simulations 
(open symbols). Closed symbols refer to the BH critical points. The case \ — 1-0 (open squares) corresponds to the Young and Adler 
molecular dynamics results^ and includes the FCC-FCC transition not considered in the patchy cases. Adapted from Refl54l 



allows to provide a prediction on the possible location of the true coexisting lines, even for coverages where MC results are not 
yet available or in regions where crystallization prevents the evaluation of the location of the metastable gas-liquid critical point. 
Note that there are no restrictions in the applicability of the BH theory, neither in thermodynamical parameters, nor in coverages, 



the only limitations being the convergence of the numerical non-linear solvers for Eqs.(44i and (45 I. On the other hand, the BH 
perturbation theory can be applied in its current form only to the one-patch case, as it lacks the possibility of discriminating the 
location of the patches. Notwithstanding these limitations, this approach is very promising as it even allows the computation of 
the fluid-solid branch, as we will see in the next section. 



E The fluid-solid coexistence 

At higher densities, a fluid-solid transition is expected in a way akin to that occurring in the fully isotropic square-well cas d 63 ! 64 ! 
For this value of the square-well amplitude A = 1.5, this transition occurs at densities that are well-separated from the fluid-fluid 
counterpart. In the isotropic SW case, this case has been studied several times with different methodology^MI In Figj9]we 
report one of these study by Young and Adler 63 (open squares) carried out using molecular dynamics techniques. In addition to 
the expected fluid-solid transition, with the solid being a face-centered-cubic (FCC) lattice, an additional FCC-FCC transition 
is visible as a large plateau at higher densities (see Fig |9|». The other open symbols in Fig|9]report results from the fluid-solid 
transition of the Kern-Frenkel model with coverages from \ = 0.9 down to the Janus limit \ — 0.5. In all cases, the final 
crystal is a FCC with 12 nearest-neighbors. There exists an additional FCC-FCC transition akin to that found in the SW case, 
corresponding to a transition from a more dilute to a denser FCC lattice, that is possible for this range of the attractive well, 
but will not be considered here. Lines in Figj9] report results from thermodynamic perturbation theory obtained by using the 
same approach outlined above for the fluid-fluid transition. In this case, a hint of the solid structural change is visible at higher 
densities, but the exact coexistence lines could not be obtained due to limitation in the convergence of the numerical algorithm 
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associated with the non-linear solver in Eqs.(44i and ( |45] l. In spite of this drawback, thermodynamic theory is able to capture 
the main features of the transition even for lower coverages (as low as \ — 0.1). It should be stressed, however, that transitions 
to different crystal structures may be envisaged at low coverages, and this feature has not be accounted in the analysis reported 
in Fig. |9] 



F Self-assembly in a predefined Kagome lattice 

The subtle interplay between the fluid-fluid and the fluid-solid transition is one of the most delicate issue in the framework of 
self-assembly patchy colloids, especially in the presence of additional effects such as inhibiting clustering transitions. This was 
already hinted before, but a very illuminating example of this is provided by the phase diagram of the triblock Janus fluid, that 
has also been considered elsewhere in this volume from the experimental point of view. 

Triblock Janus colloids are spherical colloidal particles decorated with two hydrophobic poles of tunable area, separated by 
an electrically charged middle band (triblock Janus j23. The electric charge of the particles allows for a controlled switch of 
the interaction via addition of salt, which effectively screens the overall repulsion, offering the possibility to the hydrophobic 
attraction between patches to express itself. Once deposited on a flat surface, after the addition of the salt, particles organize 
themselves into a Kagome lattice. The crystallization kinetics has been followed in real space in full details 25 . The patch width 
in the experimental system, of the order of 65 degrees corresponding to % ^ 0.57, allows for simultaneous bonding of two 



particles per patch, stabilizing the locally four-coordinated structure of the Kagome lattice (see Fig. 10 1. 

A triblock Janus particle can be modeled via the two-patches Kern-Frenkel model in which the attractive region is split in two 
parts at the opposite poles of the sphere, whereas the repulsive part is concentrated in the middle strip at the equator. The 
square-well mimics the short-range hydrophobic attraction, while the hard-sphere region represents the repulsive charge-charge 
interaction. Depending on the considered state point and on the patch amplitude several possible phases are possible, as shown 
in Ref. 27 and summarized in FigJTT] In full agreement with the experiments, we observe at comparable interaction strength the 
spontaneous nucleation of a Kagome lattice. We also find at larger pressure spontaneous crystal formation in a dense hexagonal 
structure. Such easiness to crystallize suggests that in this system crystallization barriers are comparable to the thermal energy 
at all densities. Interestingly enough, we find that for this model a (metastable) gas-liquid phase separation can be observed for 
large patch width. 

The ability of accurately describing Janus triblock particles with the KF potential is particularly rewarding and provides a strong 
support for the use of such model for predicting the self-assembly properties of this class of patchy colloids. The possibility 
of numerically exploring the sensitivity of the phase diagram to the parameters (patch width and interaction range) entering in 
the interaction potential provides an important instrument and a guide to the design of these new particles to obtain specific 
structures by self-assembly. 
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FIG. 10: From left to right: snapshot of a gas, liquid, Kagome lattice and hexagonal lattice. The Kagome and the hexagonal crystals are 
formed respectively at low and high pressures. Attractive patches are colored in blue, the hard-core remaining particle surface is colored in 
gray. Particles are free to rotate in three dimensions but are constrained to move on a flat surface. Adapted from Refl27l 
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FIG. 1 1 : Phase diagrams in the T — p for a wide (x ~ 0.57) patch model (see Ref. 27 for the short range analog). Boundaries between stable 
phases are drawn as solid black lines and metastable phase boundaries are dotted. The orange points in indicate the (metastable) gas-liquid 
critical point. The dashed line represents the metastable gas-liquid phase separation. Blue diamonds and red circles indicate the highest 
temperature at which spontaneous crystallization into the Kagome and hexagonal lattice, respectively, was detected at the corresponding 
pressure or density. Crosses indicate the coexistence points checked via direct coexistence simulations. Adapted from Refl27l 



IX Conclusions and open perspectives 

In this Chapter we have discussed some of the main theoretical approaches that have been proposed in the last few years to ad- 
dress the computation of the phase diagram in the temperature-density plane of the Kem-Frenkel model, one of the paradigmatic 
models for patchy colloids. Three different approaches, namely Monte Carlo simulations, integral equations and perturbation 
theory have been discussed, and their performances contrasted against each others. We have attempted to present the main 
ideas behind each techniques, along with some representative examples of applications. The main emphasis has been placed on 
the evaluation of the fluid-fluid (gas-liquid) and fluid-solid phase diagrams, and their importance in the framework of the self- 
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assembly processes. As a matter of fact, the exact location of the transition lines is a crucial ingredient necessary to implement 
a bottom-up engineering of materials of new generation. It is remarkable that none of these techniques are really new, being 
part of the background given in any standard graduate course in Statistical Physics. Yet, their implementation in the frame- 
work of patchy colloids often (if not always) has required significant improvements that have shed new lights on the techniques 
themselves. This is true for the newly improved Monte Carlo schemes that have been devised in this field, but also for integral 
equation and perturbation theories that have a long and glorious history. 

Several steps need still to be performed. We still need an accurate study of the role of the interaction range in selecting the 
most stable geometries. More specifically, we need to understand under which patch width and range conditions micelles and 
vesicles are the stable structures. Experiments 24 indicate that for Janus particles one dimensional structures becomes preferred 
when the interaction range is only a few percent of the particle size. We also need to develop an accurate methodology to 
predict all possible crystal structures for different types of patterned particles (including Janus), their stability fields and the 
associated nucleation rates. Understanding self-assembly into ordered pre-defined structures can indeed have relevant industrial 
applications. The Janus paradigm (and its theoretical counterpart, the Kern-Frenkel model), can play an important role in this 
process. 

The fact that notwithstanding its simplicity the Kern-Frenkel model is able to present such a complex and rich scenario in its 
phase diagram is related to the combined effect of two features. On the one hand, the short-range and reversible nature of 
the involved interactions. This allows a partial rearrangement, within a localized region in space, of the particles in search of 
the optimal minimal energy configuration, a feature that would not be possible for stronger (irreversible) covalent longer-range 
interactions. On the other hand, the specificity dictated by the patchy anisotropy is expedient in avoiding multiple degenerate 
configurations with similar energies, thus eliminating defects and polidispersity effects characteristic of isotropic colloids. 
Qualitative and semi-quantitative agreements with experiments can be obtained with the Kern-Frenkel model in some particular 
cases, and it is hoped that the contribution of this book, including both theory and experiments in a well chosen balance, will 
contribute to further strengthen this very promising route. 
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